/* 
   Trabalho 1 - Programacao Concorrente
   
   PI monte carlo versao paralela

   Henrique Teixeira Arroyo 	6419241
   Julio Polzato	        7288971
   Henry D Leão			5881632
*/

#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <pthread.h>
#include <limits.h>

#define SAMPLES 1000000000 // total samples to try (in each thread samples = SAMPLES/NUM_THREADS)
#define NUM_THREADS 50     // total threads to run		

// Global vector to store circle hits for each thread.
// This vector will be summed up later to give the total of circle hits
unsigned int in_hits[NUM_THREADS];

// myrandom implementation
const unsigned int A = 1664525, C = 1013904223;

unsigned int myrandom(unsigned int x){
        return x*A+C;
}

// function called by each thread to compute its share in the montecarlo algorithm
void *montecarlo(void *threadID){
    
    	int i,taskid, *id_ptr, iter;
    	unsigned int seed;
    	double x, y;

    	id_ptr = (int *) threadID;
    	taskid = *id_ptr;
    
	seed = myrandom(taskid*time(NULL));
	iter = SAMPLES/NUM_THREADS;
    	for(i=0;i<iter;i++){
		
		seed = myrandom(seed);
		x = (double)seed/(double)UINT_MAX;
		seed = myrandom(seed);
		y = (double)seed/(double)UINT_MAX;

		if(x*x + y*y < 1.0)
            		in_hits[taskid]++;
	}

    	pthread_exit(NULL);
}

// main function: will start the monte carlo threads and then sum up all the circle hits. With this, PI can be evaluated.
int main(){

	pthread_t run_threads;
   	unsigned int in = 0;
	double PI;
		
	int i, error, *taskids[NUM_THREADS];
    pthread_t threads[NUM_THREADS];
   
    for(i=0;i<NUM_THREADS;i++){
		taskids[i] = (int *) malloc(sizeof(int));
		*taskids[i] = i;

    	in_hits[i] = 0;
        error = pthread_create(&threads[i],NULL,montecarlo,(void *)taskids[i]);
		
		if(error){
			printf("error creating thread with id = %d \n",i);
			pthread_exit(NULL);
		}
    }	

	for(i=0;i<NUM_THREADS;i++){		    // sum up the circle hits
		pthread_join(threads[i],NULL);
		in += in_hits[i];
	}
	
	PI = (double)(4*in)/(double)SAMPLES;	    // compute PI
	
	printf("PI = %.6f\n\n",PI);
	
	return 0;
}

